#include <stdlib.h>
#include <stdio.h>
#include <string.h>
#include<suitesparse/umfpack.h>

#define MAX_LINE_SIZE 128
void sparse_add(int **J0, double **S0, int row, int col, double val);
void sparse_compress(int **J0, double **S0, int *I, int *J, double *S, int dim);

int main(int argc, char* argv[])
{
    int i,j,m,n,l;
    int *I,*J,**J0,flag;
    double val;
    double *S,**S0;
    double *rhs, *solution;

    /// this is for the smallest sparse matrix
    char* spmat_data_file = argv[1]; // "spmat.csv";
    char* rhs_data_file = argv[2]; // "rhs1.csv";
    int dim = 119425;        // 稀疏矩阵信息，这里偷懒假设事先知道！
    int nnz = 825737;

    ///
    int maxcouple = 50;   // 每行最多有几个非零元？ 给一个适当大的数，请保证数组不越界的同时内存够用!
    char dummy, strLine[MAX_LINE_SIZE];  // buffer for reading the data files
    
    /*
     * step 1. read in J0 and S0 from data file
     */
    J0 = (int**)calloc(dim, sizeof(int*));
    S0 = (double**)calloc(dim, sizeof(double*));
    for(i=0 ; i < dim; i++){
    J0[i] = (int*) calloc(maxcouple,sizeof(int));
    S0[i] = (double*) calloc(maxcouple,sizeof(double));
    for(j = 0; j < maxcouple; j++){
        J0[i][j] = 0;
        S0[i][j] = 0.0;
    }
    }

    // 打开稀疏矩阵文件 -- 假设格式知道
    FILE* fin = fopen(spmat_data_file, "r");
    if (fin != NULL)
    printf("稀疏矩阵数据文件打开成功！！\n");
    else{
    printf("稀疏矩阵数据文件打开失败，请检查！！\n");  return 0;
    }
    for (l = 0; l < nnz; l++){      // 开始读取数据
    fscanf(fin, "%d%c%d%c%lf", &i, &dummy, &j, &dummy, &val);
    //    printf("A(%d,%d) = %f\n", i-1, j-1, val);
    // add it to buffer
    sparse_add(J0, S0, i-1, j-1, val);   // please remind the -1 for index
    }
    fclose(fin);
    // for(i = 0, nnz = 0; i < dim; i++) nnz += J0[i][0]; // calculate the nnz if necesary, however, it is now assumed to be known

    /*
     * step 2. read in the rhs from dat file
    */
    rhs = (double*) calloc(dim, sizeof(double));
    fin = fopen(rhs_data_file, "r");
    if (fin != NULL)
    printf("右端项数据文件打开成功！！\n");
    else{
    printf("右端项数据文件打开失败，请检查！！\n");    return 0;
    }
    for(i = 0; i < dim; i++){
    fscanf(fin, "%le", &val);
    //      printf("rhs[%d] = %30.15f\n", i-1, val);
    rhs[i] = val;
    }
    fclose(fin);

    /*
     * step 3: (I,J,S) is the real sparse matrix data: row compressed format
     */
    I = (int*) calloc(dim+1,sizeof(int));
    J = (int*) calloc(nnz, sizeof(int));
    S = (double*)calloc(nnz, sizeof(double));
    sparse_compress(J0,S0,I,J,S,dim);

    // Finalizing the temporal data and free the memory now
    for(i = 0; i < dim; i++){
    free(J0[i]);
    free(S0[i]);
    }
    free(J0); free(S0);
    printf("Program Finished! There are %d nonzero entries, %d.\n",nnz,I[dim]);
    
    // output the sparse pattern to gnuplot data file
    /*   out = fopen("pattern.dat","w"); */
    /*   for(i = 0; i < dim; i++) */
    /*     for(l = I[i]; l < I[i+1]; l++) */
    /*       fprintf(out,"%d %d\n",J[l],-i); */
    /*   fclose(out); */

    /**
     * step 4: Solving the linear system
     */
    solution = (double*) calloc(dim, sizeof(double));
    double *null =(double*)NULL;
    void *Symbolic,*Numeric;
    (void)umfpack_di_symbolic(119425,119425,I,J,S,&Symbolic,null,null);
    (void)umfpack_di_numeric(I,J,S,Symbolic,&Numeric,null,null);
    umfpack_di_free_numeric(&Symbolic);
    umfpack_di_solve(UMFPACK_A,I,J,S,solution,rhs,Numeric,null,null);
    umfpack_di_free_numeric(&Numeric);
    free(solution);  free(rhs);
    free(I);    free(J);    free(S);

    printf("Program Finished!\n");
    return 1;

};

void sparse_add(int **J0, double **S0, int row, int col, double val)
{
  int pos;
  // find if value col is in the row's row of matrix J
  for(pos=1; pos <= J0[row][0]; pos++){
    if(J0[row][pos] == col) break;
  }
  // whether J exist or not, i is the current position to insert
  // since S is initialized to be zero, add on is safe, also the rewrite of
  // J[row][i] is also safe!! This implementation is wonderful!
  J0[row][pos] = col;
  S0[row][pos] += val;
  J0[row][0] = (pos > J0[row][0])?pos:J0[row][0];
}

void sparse_compress(int **J0, double **S0, int *I, int *J, double *S, int dim)
{
  int i,j,pos;
  // fill vector I first, whose first element must be 0
  I[0] = 0;
  for(i = 1; i <= dim; i++){
    I[i] = I[i-1] + J0[i-1][0];
  }
  // then run over J0 and S0 to fill J and S
  for(i = 0; i < dim; i++){
    for(j = 1,pos = 1; j <= J0[i][0]; j++){
      // remember that pos=0 is for the diagnal element
      if(J0[i][j] == i){
    J[I[i]] = i;
    S[I[i]] = S0[i][j];
      }
      else{
    J[I[i]+pos] = J0[i][j];
    S[I[i]+pos] = S0[i][j];
    pos++;
      }
    }
  }
}
